Fit main model (see exploratory scripts and model comparison), visualize results.
library(tidyverse)
#> Warning: package 'tidyr' was built under R version 4.0.5
library(tidylog)
library(viridis)
library(marmap)
library(rnaturalearth)
library(rnaturalearthdata)
library(sf)
#> Warning: package 'sf' was built under R version 4.0.5
library(RColorBrewer)
#> Warning: package 'RColorBrewer' was built under R version 4.0.5
library(gganimate)
library(gifski)
library(latex2exp)
library(patchwork)
library(png)
library(RCurl)
library(wesanderson)
library(qwraps2)
library(ggcorrplot)
library(ggridges)
library(brms)
#> Warning: package 'Rcpp' was built under R version 4.0.5
# remotes::install_github("pbs-assess/sdmTMB")
library(sdmTMB) # sdmTMB_0.1.4.9004
# To load entire cache in interactive r session, do:
# qwraps2::lazyload_cache_dir(path = "R/analysis/condition_model_cf_cache/html")
# Specify map ranges
ymin = 52; ymax = 60; xmin = 10; xmax = 24
map_data <- rnaturalearth::ne_countries(
scale = "medium",
returnclass = "sf", continent = "europe")
# Crop the polygon for plotting and efficiency:
# st_bbox(map_data) # find the rough coordinates
sf::sf_use_s2(FALSE)
# https://stackoverflow.com/questions/68478179/how-to-resolve-spherical-geometry-failures-when-joining-spatial-data
swe_coast <- suppressWarnings(suppressMessages(
st_crop(map_data,
c(xmin = xmin, ymin = ymin, xmax = xmax, ymax = ymax))))
# Transform our map into UTM 33 coordinates, which is the equal-area projection we fit in:
utm_zone33 <- 32633
swe_coast_proj <- sf::st_transform(swe_coast, crs = utm_zone33)
# Define plotting theme for main plot
theme_plot <- function(base_size = 11, base_family = "") {
theme_light(base_size = base_size, base_family = "") +
theme(
legend.position = "bottom",
legend.key.height = unit(0.2, "cm"),
legend.margin = margin(0, 0, 0, 0),
legend.box.margin = margin(-5, -5, -5, -5),
strip.text = element_text(size = 9, colour = 'gray10', margin = margin(b = 1, t = 1)),
strip.background = element_rect(fill = "grey95")
)
}
# Define plotting theme for facet_wrap map with years
theme_facet_map <- function(base_size = 11, base_family = "") {
theme_light(base_size = base_size, base_family = "") +
theme(
axis.text.x = element_text(angle = 90),
axis.text = element_text(size = 7),
strip.text = element_text(size = 8, colour = 'gray10', margin = margin(b = 1, t = 1)),
strip.background = element_rect(fill = "gray95"),
legend.direction = "horizontal",
legend.margin = margin(1, 1, 1, 1),
legend.box.margin = margin(0, 0, 0, 0),
legend.key.height = unit(0.4, "line"),
legend.key.width = unit(2, "line"),
legend.spacing.x = unit(0.1, 'cm'),
legend.position = "bottom",
)
}
# Make default base map plot
xmin2 <- 303000; xmax2 <- 916000; xrange <- xmax2 - xmin2
ymin2 <- 5980000; ymax2 <- 6450000; yrange <- ymax2 - ymin2
plot_map_raster <-
ggplot(swe_coast_proj) +
xlim(xmin2, xmax2) +
ylim(ymin2, ymax2) +
labs(x = "Longitude", y = "Latitude") +
geom_sf(size = 0.3) +
theme_plot() +
guides(colour = guide_colorbar(title.position = "top", title.hjust = 0.5),
fill = guide_colorbar(title.position = "top", title.hjust = 0.5))
# Read the split files and join them
d1 <- readr::read_csv("https://raw.githubusercontent.com/maxlindmark/cod-condition/master/data/for_analysis/mdat_cond_(1_2).csv")
d2 <- readr::read_csv("https://raw.githubusercontent.com/maxlindmark/cod-condition/master/data/for_analysis/mdat_cond_(2_2).csv")
d <- bind_rows(d1, d2)
unique(is.na(d$density_cod))
#> [1] FALSE
unique(is.na(d$density_cod_rec))
#> [1] FALSE TRUE
# Calculate standardized variables
d <- d %>%
mutate(ln_length_cm = log(length_cm),
ln_weight_g = log(weight_g),
year_ct = year - mean(year),
biomass_her_sc = biomass_her,
biomass_her_sd_sc = biomass_her_sd,
biomass_spr_sc = biomass_spr,
biomass_spr_sd_sc = biomass_spr_sd,
density_cod_sc = density_cod,
density_cod_rec_sc = density_cod_rec,
density_fle_sc = density_fle,
density_fle_rec_sc = density_fle_rec,
density_saduria_sc = density_saduria,
density_saduria_rec_sc = density_saduria_rec,
depth_sc = depth,
depth_rec_sc = depth_rec,
oxy_sc = oxy,
oxy_rec_sc = oxy_rec,
temp_sc = temp,
temp_rec_sc = temp_rec,
year_f = as.factor(year)) %>%
mutate_at(c("biomass_her_sc", "biomass_her_sd_sc", "biomass_spr_sc", "biomass_spr_sd_sc",
"density_cod_sc", "density_cod_rec_sc",
"density_fle_sc", "density_fle_rec_sc",
"density_saduria_sc", "density_saduria_rec_sc",
"depth_sc", "depth_rec_sc",
"oxy_sc", "oxy_rec_sc",
"temp_sc", "temp_rec_sc"
),
~(scale(.) %>% as.vector)) %>%
mutate(year = as.integer(year))
unique(is.na(d))
#> weight_g length_cm year lat lon ices_rect sub_div depth oxy temp
#> [1,] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
#> [2,] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
#> [3,] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
#> [4,] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
#> [5,] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
#> [6,] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
#> [7,] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
#> X Y density_cod density_fle density_cod_data density_fle_data
#> [1,] FALSE FALSE FALSE FALSE FALSE FALSE
#> [2,] FALSE FALSE FALSE FALSE FALSE FALSE
#> [3,] FALSE FALSE FALSE FALSE FALSE FALSE
#> [4,] FALSE FALSE FALSE FALSE TRUE TRUE
#> [5,] FALSE FALSE FALSE FALSE TRUE TRUE
#> [6,] FALSE FALSE FALSE FALSE FALSE FALSE
#> [7,] FALSE FALSE FALSE FALSE TRUE TRUE
#> density_saduria density_saduria_rec density_cod_rec density_fle_rec
#> [1,] FALSE FALSE FALSE FALSE
#> [2,] FALSE FALSE FALSE FALSE
#> [3,] FALSE FALSE FALSE FALSE
#> [4,] FALSE FALSE FALSE FALSE
#> [5,] FALSE TRUE TRUE TRUE
#> [6,] FALSE TRUE TRUE TRUE
#> [7,] FALSE FALSE FALSE FALSE
#> depth_rec oxy_rec temp_rec biomass_spr biomass_her biomass_spr_sd
#> [1,] FALSE FALSE FALSE FALSE FALSE FALSE
#> [2,] FALSE FALSE FALSE TRUE TRUE FALSE
#> [3,] FALSE FALSE FALSE TRUE TRUE TRUE
#> [4,] FALSE FALSE FALSE FALSE FALSE FALSE
#> [5,] TRUE TRUE TRUE TRUE TRUE FALSE
#> [6,] TRUE TRUE TRUE TRUE TRUE FALSE
#> [7,] FALSE FALSE FALSE TRUE TRUE FALSE
#> biomass_her_sd ln_length_cm ln_weight_g year_ct biomass_her_sc
#> [1,] FALSE FALSE FALSE FALSE FALSE
#> [2,] FALSE FALSE FALSE FALSE TRUE
#> [3,] TRUE FALSE FALSE FALSE TRUE
#> [4,] FALSE FALSE FALSE FALSE FALSE
#> [5,] FALSE FALSE FALSE FALSE TRUE
#> [6,] FALSE FALSE FALSE FALSE TRUE
#> [7,] FALSE FALSE FALSE FALSE TRUE
#> biomass_her_sd_sc biomass_spr_sc biomass_spr_sd_sc density_cod_sc
#> [1,] FALSE FALSE FALSE FALSE
#> [2,] FALSE TRUE FALSE FALSE
#> [3,] TRUE TRUE TRUE FALSE
#> [4,] FALSE FALSE FALSE FALSE
#> [5,] FALSE TRUE FALSE FALSE
#> [6,] FALSE TRUE FALSE FALSE
#> [7,] FALSE TRUE FALSE FALSE
#> density_cod_rec_sc density_fle_sc density_fle_rec_sc density_saduria_sc
#> [1,] FALSE FALSE FALSE FALSE
#> [2,] FALSE FALSE FALSE FALSE
#> [3,] FALSE FALSE FALSE FALSE
#> [4,] FALSE FALSE FALSE FALSE
#> [5,] TRUE FALSE TRUE FALSE
#> [6,] TRUE FALSE TRUE FALSE
#> [7,] FALSE FALSE FALSE FALSE
#> density_saduria_rec_sc depth_sc depth_rec_sc oxy_sc oxy_rec_sc temp_sc
#> [1,] FALSE FALSE FALSE FALSE FALSE FALSE
#> [2,] FALSE FALSE FALSE FALSE FALSE FALSE
#> [3,] FALSE FALSE FALSE FALSE FALSE FALSE
#> [4,] FALSE FALSE FALSE FALSE FALSE FALSE
#> [5,] TRUE FALSE TRUE FALSE TRUE FALSE
#> [6,] TRUE FALSE TRUE FALSE TRUE FALSE
#> [7,] FALSE FALSE FALSE FALSE FALSE FALSE
#> temp_rec_sc year_f
#> [1,] FALSE FALSE
#> [2,] FALSE FALSE
#> [3,] FALSE FALSE
#> [4,] FALSE FALSE
#> [5,] TRUE FALSE
#> [6,] TRUE FALSE
#> [7,] FALSE FALSE
unique(is.na(d$density_cod_rec))
#> [1] FALSE TRUE
unique(is.na(d$density_cod))
#> [1] FALSE
# Drop NA covariates for the correlation plot
d <- d %>% drop_na(biomass_her_sc, biomass_her_sd_sc, biomass_spr_sc, biomass_spr_sd_sc,
density_cod_sc, density_cod_rec_sc, density_fle_sc, density_fle_rec_sc,
density_saduria_sc, density_saduria_rec_sc, depth_sc, depth_rec_sc,
oxy_sc, oxy_rec_sc, temp_sc, temp_rec_sc)
# Sample size
nrow(d)
#> [1] 94295
# stdev of oxygen
sd(d$oxy)
#> [1] 1.805214
# Plot correlation between variables
d_cor <- d %>% dplyr::select("biomass_her_sc", "biomass_her_sd_sc",
"biomass_spr_sc", "biomass_spr_sd_sc",
"density_cod_sc", "density_cod_rec_sc",
"density_fle_sc", "density_fle_rec_sc",
"density_saduria_sc", "density_saduria_rec_sc",
"depth_sc", "depth_rec_sc",
"oxy_sc", "oxy_rec_sc",
"temp_sc", "temp_rec_sc")
corr <- round(cor(d_cor), 1)
ggcorrplot(corr, hc.order = TRUE, type = "lower", lab = TRUE, lab_size = 2.5) +
theme_classic() +
theme(axis.text.x = element_text(angle = 90, vjust = 0.3))
ggsave("figures/supp/condition/correlation_variables.png", width = 6.5, height = 6.5, dpi = 600)
# Plot data
d_plot <- d %>%
group_by(year, X, Y) %>%
summarise(n = n())
plot_map_raster +
geom_point(data = d_plot, aes(X*1000, Y*1000, size = n, color = n), alpha = 0.5) +
scale_color_viridis() +
scale_size(range = c(.1, 3), name = "N per haul") +
facet_wrap(~year, ncol = 5) +
guides(size = "none") +
theme_facet_map()
ggsave("figures/supp/condition/spatial_cond_data.png", width = 6.5, height = 8.5, dpi = 600)
pred_grid1 <- readr::read_csv("https://raw.githubusercontent.com/maxlindmark/cod_condition/master/data/for_analysis/pred_grid_(1_2).csv")
pred_grid2 <- readr::read_csv("https://raw.githubusercontent.com/maxlindmark/cod_condition/master/data/for_analysis/pred_grid_(2_2).csv")
pred_grid <- bind_rows(pred_grid1, pred_grid2)
unique(is.na(pred_grid$density_cod))
#> [1] FALSE
unique(is.na(pred_grid$density_cod_rec))
#> [1] FALSE
pred_grid <- pred_grid %>% mutate(year = as.integer(year),
year_f = as.factor(year))
# Scale the variables with respect to data! First calculate means in data
data_means <- d %>%
dplyr::select(biomass_her, biomass_her_sd, biomass_spr, biomass_spr_sd,
density_cod, density_cod_rec,
density_fle, density_fle_rec,
density_saduria, density_saduria_rec,
depth, depth_rec,
oxy, oxy_rec,
temp, temp_rec) %>%
mutate_all(~mean(.)) %>%
distinct(.keep_all = TRUE)
data_stdev <- d %>%
dplyr::select(biomass_her, biomass_her_sd, biomass_spr, biomass_spr_sd,
density_cod, density_cod_rec,
density_fle, density_fle_rec,
density_saduria, density_saduria_rec,
depth, depth_rec,
oxy, oxy_rec,
temp, temp_rec
) %>%
mutate_all(~sd(.)) %>%
distinct(.keep_all = TRUE)
# Before actually scaling, replace the ices-rectangle pelagic values that are NA with the mean across rectangles in the sub division.
# Replace the sub-division pelagic values that are NA with the mean in the year.
# Note that the sub-division values are not the sum of the rectangles (due to missing rectangles), so
# I need to calculate a sub-division mean across rectangles within each sub-division
pred_grid <- pred_grid %>%
group_by(year) %>%
mutate(median_biomass_sprat_across_sub_div = median(biomass_spr_sd, na.rm = TRUE),
median_biomass_herring_across_sub_div = median(biomass_her_sd, na.rm = TRUE)) %>%
ungroup() %>% # Replace sub_divsion NA's with the median across sub_divisions in that year
mutate(biomass_spr_sd = ifelse(is.na(biomass_spr_sd) == TRUE, median_biomass_sprat_across_sub_div, biomass_spr_sd),
biomass_her_sd = ifelse(is.na(biomass_her_sd) == TRUE, median_biomass_herring_across_sub_div, biomass_her_sd)) %>%
group_by(year, sub_div) %>%
mutate(median_biomass_sprat_across_rect_in_sub_div = median(biomass_spr, na.rm = TRUE),
median_biomass_herring_across_rect_in_sub_div = median(biomass_her, na.rm = TRUE)) %>%
ungroup() %>% # Replace rectangle NA's with the median across rectangles in that year and sub-division
mutate(biomass_spr = ifelse(is.na(biomass_spr) == TRUE, median_biomass_sprat_across_rect_in_sub_div, biomass_spr),
biomass_her = ifelse(is.na(biomass_her) == TRUE, median_biomass_herring_across_rect_in_sub_div, biomass_her)) %>%
# Since I still have some NAs (some sub-divisions do not have a single rectangle in some years), I will will those rectangles
# with the sub-division value divided by the number of rectangles in that sub division
group_by(year, sub_div) %>%
mutate(n_rect = length(unique(ices_rect))) %>%
ungroup() %>%
mutate(biomass_spr = ifelse(is.na(biomass_spr) == TRUE, biomass_spr_sd/n_rect, biomass_spr),
biomass_her = ifelse(is.na(biomass_her) == TRUE, biomass_her_sd/n_rect, biomass_her))
pred_grid <- pred_grid %>%
mutate(ln_length_cm = log(1),
year_ct = year - mean(year),
biomass_her_sc = (biomass_her - data_means$biomass_her) / data_stdev$biomass_her,
biomass_her_sd_sc = (biomass_her_sd - data_means$biomass_her_sd) / data_stdev$biomass_her_sd,
biomass_spr_sc = (biomass_spr - data_means$biomass_spr) / data_stdev$biomass_spr,
biomass_spr_sd_sc = (biomass_spr_sd - data_means$biomass_spr_sd) / data_stdev$biomass_spr_sd,
density_cod_sc = (density_cod - data_means$density_cod) / data_stdev$density_cod,
density_cod_rec_sc = (density_cod_rec - data_means$density_cod_rec) / data_stdev$density_cod_rec,
density_fle_sc = (density_fle - data_means$density_fle) / data_stdev$density_fle,
density_fle_rec_sc = (density_fle_rec - data_means$density_fle_rec) / data_stdev$density_fle_rec,
density_saduria_sc = (density_saduria - data_means$density_saduria) / data_stdev$density_saduria,
density_saduria_rec_sc = (density_saduria_rec - data_means$density_saduria_rec) / data_stdev$density_saduria_rec,
depth_sc = (depth - data_means$depth) / data_stdev$depth,
depth_rec_sc = (depth_rec - data_means$depth_rec) / data_stdev$depth_rec,
oxy_sc = (oxy - data_means$oxy) / data_stdev$oxy,
oxy_rec_sc = (oxy_rec - data_means$oxy_rec) / data_stdev$oxy_rec,
temp_sc = (temp - data_means$temp) / data_stdev$temp,
temp_rec_sc = (temp_rec - data_means$temp_rec) / data_stdev$temp_rec,
)
# Plot on map:
ggplot(pred_grid2, aes(X, Y)) + geom_point(size = 0.1) + facet_wrap(~year)
pred_grid %>%
drop_na() %>%
ggplot(., aes(X, Y)) + geom_point(size = 0.1) + facet_wrap(~year)
pred_grid %>%
drop_na(biomass_spr_sd) %>%
ggplot(., aes(X, Y)) + geom_point(size = 0.1) + facet_wrap(~year)
unique(is.na(pred_grid))
#> X Y depth year id oxy_q3 oxy temp_q3 temp lon lat
#> [1,] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
#> [2,] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
#> sub_div ices_rect density_saduria biomass_spr biomass_her biomass_spr_sd
#> [1,] FALSE FALSE FALSE FALSE FALSE FALSE
#> [2,] FALSE FALSE FALSE FALSE FALSE FALSE
#> biomass_her_sd density_cod density_fle depth_rec temp_rec oxy_rec
#> [1,] FALSE FALSE FALSE FALSE FALSE FALSE
#> [2,] FALSE FALSE FALSE FALSE FALSE FALSE
#> density_cod_rec density_fle_rec density_saduria_rec year_f
#> [1,] FALSE FALSE FALSE FALSE
#> [2,] FALSE FALSE FALSE FALSE
#> median_biomass_sprat_across_sub_div median_biomass_herring_across_sub_div
#> [1,] FALSE FALSE
#> [2,] FALSE FALSE
#> median_biomass_sprat_across_rect_in_sub_div
#> [1,] FALSE
#> [2,] TRUE
#> median_biomass_herring_across_rect_in_sub_div n_rect ln_length_cm year_ct
#> [1,] FALSE FALSE FALSE FALSE
#> [2,] TRUE FALSE FALSE FALSE
#> biomass_her_sc biomass_her_sd_sc biomass_spr_sc biomass_spr_sd_sc
#> [1,] FALSE FALSE FALSE FALSE
#> [2,] FALSE FALSE FALSE FALSE
#> density_cod_sc density_cod_rec_sc density_fle_sc density_fle_rec_sc
#> [1,] FALSE FALSE FALSE FALSE
#> [2,] FALSE FALSE FALSE FALSE
#> density_saduria_sc density_saduria_rec_sc depth_sc depth_rec_sc oxy_sc
#> [1,] FALSE FALSE FALSE FALSE FALSE
#> [2,] FALSE FALSE FALSE FALSE FALSE
#> oxy_rec_sc temp_sc temp_rec_sc
#> [1,] FALSE FALSE FALSE
#> [2,] FALSE FALSE FALSE
pred_grid %>%
drop_na(biomass_spr_sc, biomass_spr_sd_sc, biomass_her_sc, biomass_her_sd_sc) %>%
ggplot(., aes(X, Y)) + geom_point(size = 0.1) + facet_wrap(~year)
spde <- make_mesh(d, xy_cols = c("X", "Y"),
n_knots = 100,
type = "kmeans", seed = 42)
# Plot and save spde
png(file = "figures/supp/condition/spde.png", units = "in", width = 6.5, height = 6.5, res = 300)
plot(spde)
dev.off()
ptm <- proc.time()
mnull <- sdmTMB(
formula = ln_weight_g ~ ln_length_cm,
data = d,
time = NULL,
mesh = spde,
family = student(link = "identity", df = 5),
spatiotemporal = "off",
spatial = "off",
silent = TRUE,
reml = FALSE
)
proc.time() - ptm
#> user system elapsed
#> 1.611 0.116 1.763
sanity(mnull)
#> ✖ Non-linear minimizer did not converge: do not trust this model!
#> ℹ Try simplifying the model, adjusting the mesh, or adding priors
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigenvalues detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✔ No standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
# mlm <- lm(
# formula = ln_weight_g ~ ln_length_cm,
# data = d)
# bm <- brms::brm(formula = ln_weight_g ~ ln_length_cm, data = d, family = student())
# summary(bm)
# plot(bm)
# brms::pp_check(bm)
# Right, so even though sdmTMB warns, the actual estimates are very similar to lm and brms, student-t or gaussian. So will just play along with it
# Extract average a and b across the total dataset
a <- as.numeric(tidy(mnull)[1, "estimate"])
b <- as.numeric(tidy(mnull)[2, "estimate"])
a
#> [1] -4.638609
b
#> [1] 2.984044
d$weight_g_avg_pred <- exp(a)*d$length_cm^b
plot(d$weight_g ~ d$length_cm)
plot(d$weight_g_avg_pred ~ d$length_cm)
plot(d$weight_g_avg_pred ~ d$weight_g); abline(a = 0, b = 1, col = "red")
d$le_cren <- d$weight_g / d$weight_g_avg_pred
d$log_le_cren <- log(d$le_cren)
hist(d$le_cren)
ptm <- proc.time()
mfull <- sdmTMB(
formula = log_le_cren ~ -1 + year_f + biomass_her_sc + biomass_her_sd_sc +
biomass_spr_sc + biomass_spr_sd_sc +
density_cod_sc + density_cod_rec_sc +
density_fle_sc + density_fle_rec_sc +
density_saduria_sc + density_saduria_rec_sc +
depth_sc + depth_rec_sc +
oxy_sc + oxy_rec_sc +
temp_sc + temp_rec_sc,
data = d,
time = "year",
mesh = spde,
family = student(link = "identity", df = 5),
spatiotemporal = "iid",
spatial = "on",
silent = TRUE,
reml = FALSE,
control = sdmTMBcontrol(newton_loops = 1)#,
# priors = sdmTMBpriors(b = normal(c(rep(NA, 27), rep(0, 16)),
# c(rep(NA, 27), rep(1, 16))))
)
system("say Just finished!")
proc.time() - ptm
#> user system elapsed
#> 883.707 43.842 944.490
sanity(mfull)
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigenvalues detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✔ No standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
tidy(mfull, conf.int = TRUE) %>% arrange(desc(abs(estimate)))
tidy(mfull, effects = "ran_par", conf.int = TRUE)
# Arrange coefficients
tidy(mfull, conf.int = TRUE) %>% filter(!grepl("year", term)) %>% arrange(desc(estimate))
#> filter: removed 27 rows (63%), 16 rows remaining
tidy(mfull, conf.int = TRUE) %>% filter(!grepl("year", term)) %>% arrange(estimate)
#> filter: removed 27 rows (63%), 16 rows remaining
# Average magnitude of fixed effects?
tidy(mfull, conf.int = TRUE) %>%
filter(!grepl("year", term)) %>%
mutate(abs_coef = abs(estimate)) %>%
summarise(mean_abs_coef = mean(abs_coef))
#> filter: removed 27 rows (63%), 16 rows remaining
#> mutate: new variable 'abs_coef' (double) with 16 unique values and 0% NA
#> summarise: now one row and one column, ungrouped
tidy(mfull, conf.int = TRUE, effects = "ran_par") %>%
filter(term %in% c("sigma_O", "sigma_E")) %>%
summarise(mean_estimate = mean(estimate))
#> filter: removed 2 rows (50%), 2 rows remaining
#> summarise: now one row and one column, ungrouped
# This model is for the sensitivity analysis, comparing estimates between different model subsets etc. These tend to fit better with only a spatiotemporal term (spatial stdev tends to blow up), so I want them to be as comparable as possible
ptm <- proc.time()
mfull_spatiotemporal <- sdmTMB(
formula = log_le_cren ~ -1 + year_f + biomass_her_sc + biomass_her_sd_sc +
biomass_spr_sc + biomass_spr_sd_sc +
density_cod_sc + density_cod_rec_sc +
density_fle_sc + density_fle_rec_sc +
density_saduria_sc + density_saduria_rec_sc +
depth_sc + depth_rec_sc +
oxy_sc + oxy_rec_sc +
temp_sc + temp_rec_sc,
data = d,
time = "year",
mesh = spde,
family = student(link = "identity", df = 5),
spatiotemporal = "iid",
spatial = "off",
silent = TRUE,
reml = FALSE,
control = sdmTMBcontrol(newton_loops = 1)
)
system("say Just finished!")
proc.time() - ptm
#> user system elapsed
#> 672.519 20.072 704.688
sanity(mfull_spatiotemporal)
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigenvalues detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✔ No standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
tidy(mfull_spatiotemporal, conf.int = TRUE) %>% arrange(desc(abs(estimate)))
tidy(mfull_spatiotemporal, effects = "ran_par", conf.int = TRUE)
# Extract random and fixed coefficients from the full model
mfull_est_st <- bind_rows(tidy(mfull_spatiotemporal, effects = "ran_par", conf.int = TRUE) %>%
filter(term == "sigma_E"),
tidy(mfull_spatiotemporal, effects = "fixed", conf.int = TRUE) %>%
filter(!grepl('year', term))) %>%
mutate(term = factor(term))
#> filter: removed 2 rows (67%), one row remaining
#> filter: removed 27 rows (63%), 16 rows remaining
#> mutate: converted 'term' from character to factor (0 new NA)
# Extract coefficients and save as csv
mfull_est_st <- mfull_est_st %>%
mutate(group = "Herring",
group = ifelse(term %in% c("biomass_spr_sc", "biomass_spr_sd_sc"), "Sprat", group),
group = ifelse(term %in% c("density_cod_rec_sc", "density_cod_sc"), "Cod", group),
group = ifelse(term %in% c("density_fle_rec_sc", "density_fle_sc"), "Flounder", group),
group = ifelse(term %in% c("density_saduria_rec_sc", "density_saduria_sc"), "Saduria", group),
group = ifelse(term %in% c("depth_rec_sc", "depth_sc"), "Depth", group),
group = ifelse(term %in% c("oxy_rec_sc", "oxy_sc"), "Oxygen", group),
group = ifelse(term %in% c("temp_rec_sc", "temp_sc"), "Temp", group),
group = ifelse(term == "sigma_E", "Random", group),
) %>%
mutate(scale = "Subdivision",
scale = ifelse(term %in% c("biomass_her_sc", "biomass_spr_sc", "density_cod_rec_sc", "density_fle_rec_sc",
"density_saduria_rec_sc", "depth_rec_sc", "oxy_rec_sc", "temp_rec_sc"), "Rectangle", scale),
scale = ifelse(term %in% c("density_cod_sc", "density_fle_sc", "density_saduria_sc",
"depth_sc", "oxy_sc", "temp_sc"), "Haul", scale),
scale = ifelse(term == "sigma_E", "Spatial/spatiotemporal s.d.", scale),
) %>%
mutate(term2 = ifelse(term == "sigma_E", 2, 1))
#> mutate: new variable 'group' (character) with 9 unique values and 0% NA
#> mutate: new variable 'scale' (character) with 4 unique values and 0% NA
#> mutate: new variable 'term2' (double) with 2 unique values and 0% NA
# Now save this for comparison in sensitivity analysis script
write.csv(mfull_est_st, "output/mfull_est_st.csv")
# Save model (so that I can add density predictions to the pred grid for condition)
saveRDS(mfull, "output/mfull.rds")
mcmc_res <- residuals(mfull, type = "mle-mcmc", mcmc_iter = 5000, mcmc_warmup = 2000)
png(file = "figures/supp/condition/qq.png", units = "in", width = 6.5, height = 6.5, res = 300)
qqnorm(mcmc_res);qqline(mcmc_res)
dev.off()
d$residuals <- mcmc_res[, 1]
r2.sdmTMBb_fe <- tidy(mfull)
b_fe <- stats::setNames(b_fe$estimate, b_fe$term)
X <- mfull$tmb_data$X_ij
X <- matrix(unlist(X), ncol = length(b_fe))
VarF <- var(as.vector(b_fe %*% t(X))) # variance from fixed-effects
# What if I skip the year effect?
b_fe2 <- tidy(mfull) %>% filter(!grepl('year', term))
#> filter: removed 27 rows (63%), 16 rows remaining
b_fe2 <- stats::setNames(b_fe2$estimate, b_fe2$term)
X2 <- mfull$tmb_data$X_ij
first_index <- match('year_f2019', names(as.data.frame(mfull$tmb_data$X_ij[[1]]))) + 1
last_index <- ncol(mfull$tmb_data$X_ij[[1]])
X2 <- matrix(unlist(X2), ncol = length(b_fe))[, first_index:last_index] # skip year columns
VarF2 <- var(as.vector(b_fe2 %*% t(X2))) # variance from fixed-effects
b_ran <- tidy(mfull, "ran_par")
sigma_O <- b_ran$estimate[b_ran$term == "sigma_O"] # spatial standard deviation
sigma_E <- b_ran$estimate[b_ran$term == "sigma_E"] # spatiotemporal standard deviation
VarO <- sigma_O^2 # spatial variance
VarE <- sigma_E^2 # spatiotemporal variance
# Calculate obs - pred
d_r2 <- d
d_r2$pred <- predict(mfull)
d_r2$resid_manual <- d_r2$le_cren - d_r2$pred$est
#hist(d_r2$resid_manual)
VarR <- var(as.vector(d_r2$resid_manual)) # "residual" variance
# Marginal R2s
denominator <- VarF + VarO + VarE + VarR
denominator_no_yr <- VarF2 + VarO + VarE + VarR
marg <- VarF/denominator
marg_no_yr <- VarF2/denominator_no_yr
out <- data.frame(
marginal = marg,
partial_spatial = VarO/denominator,
partial_spatiotemporal = VarE/denominator,
conditional = (VarF + VarO + VarE)/denominator,
all_random = (VarO + VarE)/denominator,
marginal_random_ratio = marg / ((VarO + VarE)/denominator),
random_marginal_ratio = ((VarO + VarE)/denominator) / marg
)
out
out_no_yr <- data.frame(
marginal = marg_no_yr,
partial_spatial = VarO/denominator_no_yr,
partial_spatiotemporal = VarE/denominator_no_yr,
conditional = (VarF2 + VarO + VarE)/denominator_no_yr,
all_random = (VarO + VarE)/denominator_no_yr,
marginal_random_ratio = marg_no_yr / ((VarO + VarE)/denominator_no_yr),
random_marginal_ratio = ((VarO + VarE)/denominator_no_yr) / marg_no_yr
)
out_no_yr
# Residuals on map
plot_map_raster +
geom_point(data = d, aes(x = X * 1000, y = Y * 1000, color = residuals), size = 0.5) +
scale_colour_gradient2() +
facet_wrap(~ year, ncol = 5) +
labs(color = "residuals") +
theme_facet_map() +
geom_sf(size = 0.1)
ggsave("figures/supp/condition/spatial_resid.png", width = 6.5, height = 8.5, dpi = 600)
# Residuals vs length and year
ggplot(d, aes(x = length_cm, y = residuals, color = residuals)) +
geom_point(size = 0.1, alpha = 1) +
geom_hline(yintercept = 0, size = 0.5, color = "black", linetype = 2, alpha = 0.5) +
labs(x = "length [cm]", y = "Residuals", color = "") +
stat_smooth(color = "black", size = 0.5) +
facet_wrap(~year, ncol = 5) +
scale_color_gradient2() +
theme_plot()
ggsave("figures/supp/condition/residuals_vs_length_and_year.png", width = 6.5, height = 8.5, dpi = 600)
ggplot(d, aes(year, length_cm, z = residuals)) +
stat_summary_2d(bins = 40) +
scale_fill_gradient2()
# ggsave("figures/supp/condition/residuals_vs_length_and_year2.png", width = 6.5, height = 6.5, dpi = 600)
get_index_simspred_avg_sim <- predict(mfull, newdata = pred_grid %>% filter(depth > 10 & depth < 110), nsim = 500)
#> filter: removed 50,976 rows (20%), 207,603 rows remaining
#ncells <- filter(pred_grid, year == max(pred_grid$year)) %>% nrow()
# index_avg_sim <- get_index_sims(pred_avg_sim, area = rep(1/ncells, nrow(pred_avg_sim)))
# index_avg_sims <- get_index_sims(pred_avg_sim, area = rep(1/ncells, nrow(pred_avg_sim)),
# return_sims = TRUE) # This is for calculation on samples!
# Old version above: now we do a weighted version. We could do it directly in the get_index_sims or manually in the prediction grid (but then we don't get an CIs)
# This is so that we can divide by the sum of weights (needs to be done by year)
weight_sum <- pred_grid %>%
filter(depth > 10 & depth < 110) %>%
group_by(year) %>%
summarise(density_cod = sum(density_cod))
#> filter: removed 50,976 rows (20%), 207,603 rows remaining
#> group_by: one grouping variable (year)
#> summarise: now 27 rows and 2 columns, ungrouped
index_avg_sim <- get_index_sims(pred_avg_sim,
area = filter(pred_grid, depth > 10 & depth < 110)$density_cod) %>%
left_join(weight_sum) %>%
mutate(
est = est / density_cod,
lwr = lwr / density_cod,
upr = upr / density_cod
)
#> filter: removed 50,976 rows (20%), 207,603 rows remaining
#> We generally recommend using `get_index(..., bias_correct = TRUE)`
#> rather than `get_index_sims()`.
#> Warning: Default `agg_function` and `area_function` apply to
#> values provided that are in log space. This matrix may be of
#> `type` = `response` or in a different link space.
#> Joining, by = "year"
#> left_join: added one column (density_cod)
#> > rows only in x 0
#> > rows only in y ( 0)
#> > matched rows 27
#> > ====
#> > rows total 27
#> mutate: changed 27 values (100%) of 'est' (0 new NA)
#> changed 27 values (100%) of 'lwr' (0 new NA)
#> changed 27 values (100%) of 'upr' (0 new NA)
index_avg_sims <- get_index_sims(pred_avg_sim,
area = filter(pred_grid, depth > 10 & depth < 110)$density_cod,
return_sims = TRUE) %>%
left_join(weight_sum) %>%
mutate(est = .value / density_cod)
#> filter: removed 50,976 rows (20%), 207,603 rows remaining
#> We generally recommend using `get_index(..., bias_correct = TRUE)`
#> rather than `get_index_sims()`.
#> Warning: Default `agg_function` and `area_function` apply to
#> values provided that are in log space. This matrix may be of
#> `type` = `response` or in a different link space.
#> Joining, by = "year"
#> left_join: added one column (density_cod)
#> > rows only in x 0
#> > rows only in y ( 0)
#> > matched rows 13,500
#> > ========
#> > rows total 13,500
#> mutate: new variable 'est' (double) with 13,500 unique values and 0% NA
ggplot(index_avg_sim, aes(year, est, ymin = lwr, ymax = upr)) +
geom_ribbon(alpha = 0.2, color = NA) +
geom_line(size = 1) +
geom_line(data = filter(index_avg_sims, .iteration < 26),
aes(year, est, group = .iteration), inherit.aes = FALSE, alpha = 0.3)
#> filter: removed 12,825 rows (95%), 675 rows remaining
pred_grid_fe <- pred_grid %>%
mutate(year_f = as.factor(1993))
#> mutate: changed 249,002 values (96%) of 'year_f' (0 new NA)
pred_avg_sim_fe <- predict(mfull, newdata = pred_grid_fe %>% filter(depth > 10 & depth < 110), nsim = 500, re_form = NA)
#> filter: removed 50,976 rows (20%), 207,603 rows remaining
weight_sum <- pred_grid %>%
filter(depth > 10 & depth < 110) %>%
group_by(year) %>%
summarise(density_cod = sum(density_cod))
#> filter: removed 50,976 rows (20%), 207,603 rows remaining
#> group_by: one grouping variable (year)
#> summarise: now 27 rows and 2 columns, ungrouped
index_avg_sim_fe <- get_index_sims(pred_avg_sim_fe,
area = filter(pred_grid, depth > 10 & depth < 110)$density_cod) %>%
left_join(weight_sum) %>%
mutate(
est = est / density_cod,
lwr = lwr / density_cod,
upr = upr / density_cod
)
#> filter: removed 50,976 rows (20%), 207,603 rows remaining
#> We generally recommend using `get_index(..., bias_correct = TRUE)`
#> rather than `get_index_sims()`.
#> Warning: Default `agg_function` and `area_function` apply to
#> values provided that are in log space. This matrix may be of
#> `type` = `response` or in a different link space.
#> Joining, by = "year"
#> left_join: added one column (density_cod)
#> > rows only in x 0
#> > rows only in y ( 0)
#> > matched rows 27
#> > ====
#> > rows total 27
#> mutate: changed 27 values (100%) of 'est' (0 new NA)
#> changed 27 values (100%) of 'lwr' (0 new NA)
#> changed 27 values (100%) of 'upr' (0 new NA)
index_avg_sims_fe <- get_index_sims(pred_avg_sim_fe,
area = filter(pred_grid, depth > 10 & depth < 110)$density_cod,
return_sims = TRUE) %>%
left_join(weight_sum) %>%
mutate(est = .value / density_cod)
#> filter: removed 50,976 rows (20%), 207,603 rows remaining
#> We generally recommend using `get_index(..., bias_correct = TRUE)`
#> rather than `get_index_sims()`.
#> Warning: Default `agg_function` and `area_function` apply to
#> values provided that are in log space. This matrix may be of
#> `type` = `response` or in a different link space.
#> Joining, by = "year"
#> left_join: added one column (density_cod)
#> > rows only in x 0
#> > rows only in y ( 0)
#> > matched rows 13,500
#> > ========
#> > rows total 13,500
#> mutate: new variable 'est' (double) with 13,500 unique values and 0% NA
# With random effects but no fixed year effect
# pred_grid_re <- pred_grid %>%
# mutate(year_f = as.factor(1993))
#
# pred_avg_sim_re <- predict(mfull, newdata = pred_grid_re %>% filter(depth > 10 & depth < 110), nsim = 500)
#
# weight_sum <- pred_grid %>%
# filter(depth > 10 & depth < 110) %>%
# group_by(year) %>%
# summarise(density_cod = sum(density_cod))
#
# index_avg_sim_re <- get_index_sims(pred_avg_sim_re,
# area = filter(pred_grid, depth > 10 & depth < 110)$density_cod) %>%
# left_join(weight_sum) %>%
# mutate(
# est = est / density_cod,
# lwr = lwr / density_cod,
# upr = upr / density_cod
# )
#
# index_avg_sims_re <- get_index_sims(pred_avg_sim_re,
# area = filter(pred_grid, depth > 10 & depth < 110)$density_cod,
# return_sims = TRUE) %>%
# left_join(weight_sum) %>%
# mutate(est = .value / density_cod)
# Plot all together
index_comp_sim <- bind_rows(#index_avg_sim_re %>% mutate(Prediction = "Random & Fixed effects (no yr)"),
index_avg_sim_fe %>% mutate(Prediction = "Fixed effects (no yr)"),
index_avg_sim %>% mutate(Prediction = "Full"))
#> mutate: new variable 'Prediction' (character) with one unique value and 0% NA
#> mutate: new variable 'Prediction' (character) with one unique value and 0% NA
index_comp_sims <- bind_rows(#index_avg_sims_re %>% mutate(Prediction = "Random & Fixed effects (no yr)"),
index_avg_sims_fe %>% mutate(Prediction = "Fixed effects (no yr)"),
index_avg_sims %>% mutate(Prediction = "Full")) %>%
filter(.iteration < 26)
#> mutate: new variable 'Prediction' (character) with one unique value and 0% NA
#> mutate: new variable 'Prediction' (character) with one unique value and 0% NA
#> filter: removed 25,650 rows (95%), 1,350 rows remaining
ggplot(index_comp_sim, aes(year, est, ymin = lwr, ymax = upr, color = Prediction, fill = Prediction)) +
geom_ribbon(alpha = 0.2, color = NA) +
geom_line(size = 1) +
geom_line(data = filter(index_comp_sims, .iteration < 26),
aes(year, est, group = interaction(.iteration, Prediction), color = Prediction), inherit.aes = FALSE, alpha = 0.3) +
NULL
#> filter: no rows removed
# Reshape data again to calculate row-wise differences
year_diff <- index_avg_sims %>%
dplyr::select(-.value, -density_cod) %>%
filter(year %in% c(1993, 2019)) %>%
pivot_wider(names_from = year, values_from = est) %>%
mutate("2019-1993" = `2019` - `1993`) %>%
pivot_longer(2:4) %>%
ungroup()
#> filter: removed 12,500 rows (93%), 1,000 rows remaining
#> pivot_wider: reorganized (year, est) into (1993, 2019) [was 1000x3, now 500x3]
#> mutate: new variable '2019-1993' (double) with 500 unique values and 0% NA
#> pivot_longer: reorganized (1993, 2019, 2019-1993) into (name, value) [was 500x4, now 1500x3]
#> ungroup: no grouping variables
# Plot these
p1 <- ggplot(year_diff) +
geom_density(data = filter(year_diff, !name == "2019-1993"),
aes(value, fill = name), alpha = 0.6) +
scale_fill_brewer(palette = "Dark2") +
coord_cartesian(expand = 0)
#> filter: removed 500 rows (33%), 1,000 rows remaining
p2 <- ggplot(year_diff) +
geom_density(data = filter(year_diff, name == "2019-1993"),
aes(value), alpha = 0.6, fill = "grey") +
geom_vline(xintercept = 0, linetype = 2)
#> filter: removed 1,000 rows (67%), 500 rows remaining
p1 / p2
# Calculate quantiles for each level (to go in the results section)
year_diff %>%
group_by(name) %>%
summarise(median = quantile(value, probs = 0.5),
upr97.5 = quantile(value, probs = 0.975),
lwr2.5 = quantile(value, probs = 0.025))
#> group_by: one grouping variable (name)
#> summarise: now 3 rows and 4 columns, ungrouped
# Percent change in condition factor
year_diff %>%
group_by(name) %>%
pivot_wider(names_from = name, values_from = value) %>%
mutate(perc_change = ((`2019`-`1993`)/`1993`)*100) %>%
ungroup() %>%
summarise(lwr2.5 = quantile(perc_change, probs = 0.025),
median = quantile(perc_change, probs = 0.5),
upr97.5 = quantile(perc_change, probs = 0.975))
#> group_by: one grouping variable (name)
#> pivot_wider: reorganized (name, value) into (1993, 2019, 2019-1993) [was 1500x3, now 500x4]
#> mutate: new variable 'perc_change' (double) with 500 unique values and 0% NA
#> ungroup: no grouping variables
#> summarise: now one row and 3 columns, ungrouped
Evaluate sensitivity to using a pred grid with a subset excluding low density areas
# Now do the same but excluding the areas not sampled in the pred grid
# ncells_sub <- filter(pred_grid, year == max(pred_grid$year) & depth > 40 & depth < 120) %>% nrow()
# pred_grid_sub <- filter(pred_grid, depth > 40 & depth < 120)
#
# pred_avg_sim_sub <- predict(mfull, newdata = pred_grid_sub, nsim = 100)
#
# index_avg_sim_sub <- get_index_sims(pred_avg_sim_sub, area = rep(1/ncells_sub, nrow(pred_avg_sim_sub)))
#
# # Combine and plot & compare
# index_avg_sim_comp <- bind_rows(mutate(index_avg_sim, area = "full"),
# mutate(index_avg_sim_sub, area = "subset"))
#
# ggplot(index_avg_sim_comp, aes(year, est, ymin = lwr, ymax = upr, color = area, fill = area)) +
# geom_line() +
# geom_ribbon(alpha = 0.4, color = NA)
div_index_list <- list()
for(i in unique(pred_grid$sub_div)){
pred_grid_sub <- pred_grid %>% filter(sub_div == i & depth > 10 & depth < 110)
pred_avg_sim_sub <- predict(mfull, newdata = pred_grid_sub, nsim = 500)
weight_sum_sub <- pred_grid_sub %>%
group_by(year) %>%
summarise(density_cod = sum(density_cod))
index_avg_sim <- get_index_sims(pred_avg_sim_sub,
area = pred_grid_sub$density_cod) %>%
left_join(weight_sum_sub) %>%
mutate(est = est / density_cod,
lwr = lwr / density_cod,
upr = upr / density_cod) %>%
mutate(sub_div = i)
div_index_list[[i]] <- index_avg_sim
}
#> Warning: Default `agg_function` and `area_function` apply to
#> values provided that are in log space. This matrix may be of
#> `type` = `response` or in a different link space.
#> Default `agg_function` and `area_function` apply to
#> values provided that are in log space. This matrix may be of
#> `type` = `response` or in a different link space.
#> Default `agg_function` and `area_function` apply to
#> values provided that are in log space. This matrix may be of
#> `type` = `response` or in a different link space.
#> Default `agg_function` and `area_function` apply to
#> values provided that are in log space. This matrix may be of
#> `type` = `response` or in a different link space.
#> Default `agg_function` and `area_function` apply to
#> values provided that are in log space. This matrix may be of
#> `type` = `response` or in a different link space.
div_index <- bind_rows(div_index_list)
# Spatiotemporally standardized indicies
div_index2 <- bind_rows(div_index %>% dplyr::select(year, est, lwr, upr, sub_div) %>% mutate(sub_div = as.factor(sub_div)),
mutate(index_avg_sim, sub_div = as.factor("Full model prediction")))
# Add in raw mean of data:
div_index2 <- bind_rows(div_index2,
d %>%
group_by(year) %>%
summarise(est = exp(mean(log_le_cren))) %>%
mutate(sub_div = "Empirical mean"))
# Add fixed effect predictions
div_index2 <- bind_rows(div_index2,
index_avg_sim_fe %>% mutate(sub_div = "Fixed effects (year omitted)"))
# Add sims for total and fixed effect prediction
index_comp_sims <- bind_rows(index_avg_sims_fe %>% mutate(sub_div = "Fixed effects (year omitted)"),
index_avg_sims %>% mutate(sub_div = "Full model prediction")) %>%
filter(.iteration < 26)
pal <- brewer.pal(n = 4, name = "Set1")[2:4]
# Plot!
p1 <- div_index2 %>%
filter(sub_div %in% c("Full model prediction", "Empirical mean", "Fixed effects (year omitted)")) %>%
mutate(plot_facet = "Total and data") %>%
mutate(plot_group = ifelse(sub_div == "Empirical mean", 1, 0)) %>%
ggplot(aes(year, est, color = reorder(sub_div, plot_group))) +
labs(y = "Le Cren's condition factor", x = "", color = "") +
geom_ribbon(aes(x = year, y = est, ymin = lwr, ymax = upr, fill = sub_div), color = NA, alpha = 0.25) +
geom_line(data = filter(index_comp_sims, .iteration < 26),
aes(year, est, group = interaction(.iteration, sub_div), color = sub_div), inherit.aes = FALSE, alpha = 0.3) +
geom_line(size = 0.8) +
facet_wrap(~plot_facet, ncol = 1) +
scale_color_brewer(palette = "Set2") +
scale_fill_brewer(palette = "Set2") +
scale_x_continuous(breaks = c(1994, 2000, 2006, 2012, 2018)) +
# scale_fill_manual(values = pal) +
# scale_color_manual(values = pal) +
theme_plot(base_size = 10) +
theme(plot.margin = unit(c(0.4, 0.4, 0, 0.4), "cm"),
legend.position = c(0.22, 0.1),
legend.direction = "horizontal",
legend.spacing.y = unit(0, 'cm'),
legend.key.size = unit(1, "cm"),
legend.key.width = unit(0.5, "cm"),
legend.title = element_text(size = 8),
legend.text = element_text(size = 7),
legend.background = element_rect(fill = NA),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank()) +
guides(fill = "none",
color = guide_legend(ncol = 1, title.position = "top", title.hjust = 0.5))
p1
#> Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning
#> -Inf
p2 <- div_index2 %>%
filter(!sub_div %in% c("Full model prediction", "Empirical mean", "Fixed effects (year omitted)")) %>%
mutate(plot_facet = "Subdivision") %>%
ggplot(aes(year, est, color = sub_div, fill = sub_div)) +
scale_x_continuous(breaks = scales::pretty_breaks(n = 10)) +
labs(y = "Le Cren's condition factor", x = "Year") +
#geom_ribbon(aes(x = year, y = est, ymin = lwr, ymax = upr), color = NA, alpha = 0.15) +
geom_line() +
facet_wrap(~plot_facet, ncol = 1) +
scale_color_brewer(palette = "Dark2", name = "") +
scale_fill_brewer(palette = "Dark2") +
scale_x_continuous(breaks = c(1994, 2000, 2006, 2012, 2018)) +
theme_plot(base_size = 10) +
theme(plot.margin = unit(c(0, 0.4, 0.4, 0.4), "cm"),
legend.position = c(0.92, 0.88),
legend.direction = "horizontal",
legend.spacing.y = unit(0, 'cm'),
legend.key.size = unit(0.8, "cm"),
legend.key.width = unit(0.3, "cm"),
legend.title = element_text(size = 8),
legend.text = element_text(size = 7),
legend.background = element_rect(fill = NA),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank()) +
guides(color = guide_legend(ncol = 1, title.position = "top", title.hjust = 0.5),
fill = "none")
p2
p1 / p2 + plot_annotation(tag_levels = "A")
#> Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning
#> -Inf
#ggsave("figures/cond_index.pdf", width = 11, height = 17, units = "cm")
# What is one standard deviation of oxygen?
sd(d$oxy)
#> [1] 1.805214
# Extract random and fixed coefficients from the full model
mfull_est <- bind_rows(tidy(mfull, effects = "ran_par", conf.int = TRUE) %>%
filter(term %in% c("sigma_O", "sigma_E")),
tidy(mfull, effects = "fixed", conf.int = TRUE) %>%
filter(!grepl('year', term))) %>%
mutate(term = factor(term))
mfull_est <- mfull_est %>%
mutate(group = "Herring",
group = ifelse(term %in% c("biomass_spr_sc", "biomass_spr_sd_sc"), "Sprat", group),
group = ifelse(term %in% c("density_cod_rec_sc", "density_cod_sc"), "Cod", group),
group = ifelse(term %in% c("density_fle_rec_sc", "density_fle_sc"), "Flounder", group),
group = ifelse(term %in% c("density_saduria_rec_sc", "density_saduria_sc"), "Saduria", group),
group = ifelse(term %in% c("depth_rec_sc", "depth_sc"), "Depth", group),
group = ifelse(term %in% c("oxy_rec_sc", "oxy_sc"), "Oxygen", group),
group = ifelse(term %in% c("temp_rec_sc", "temp_sc"), "Temp", group),
group = ifelse(term == "sigma_E", "Random", group),
group = ifelse(term == "sigma_O", "Random", group)
)
mfull_est <- mfull_est %>%
mutate(scale = "Subdivision",
scale = ifelse(term %in% c("biomass_her_sc", "biomass_spr_sc", "density_cod_rec_sc", "density_fle_rec_sc",
"density_saduria_rec_sc", "depth_rec_sc", "oxy_rec_sc", "temp_rec_sc"), "Rectangle", scale),
scale = ifelse(term %in% c("density_cod_sc", "density_fle_sc", "density_saduria_sc",
"depth_sc", "oxy_sc", "temp_sc"), "Haul", scale),
scale = ifelse(term == "sigma_E", "Spatial/spatiotemporal s.d.", scale),
scale = ifelse(term == "sigma_O", "Spatial/spatiotemporal s.d.", scale)
)
# Quick calculation:
mfull_est %>%
mutate(par_type = ifelse(term %in% c("sigma_O", "sigma_E"), "re", "fe")) %>%
group_by(par_type) %>%
summarise(mean_est = mean(estimate)) %>%
pivot_wider(names_from = par_type, values_from = mean_est) %>%
mutate(ratio = re/fe)
# Plot effects
# This is the order changed labels should go in!
levels(mfull_est$term)
#> [1] "biomass_her_sc" "biomass_her_sd_sc" "biomass_spr_sc"
#> [4] "biomass_spr_sd_sc" "density_cod_rec_sc" "density_cod_sc"
#> [7] "density_fle_rec_sc" "density_fle_sc" "density_saduria_rec_sc"
#> [10] "density_saduria_sc" "depth_rec_sc" "depth_sc"
#> [13] "oxy_rec_sc" "oxy_sc" "sigma_E"
#> [16] "sigma_O" "temp_rec_sc" "temp_sc"
# I want the random effects to be dark gray...
pal <- brewer.pal(n = 8, name = "Dark2")
pal2 <- c(pal[1:5], "black", pal[6:8])
# Sort the terms so that random effects are at the top...
mfull_est <- mfull_est %>%
mutate(term2 = ifelse(term %in% c("sigma_E", "sigma_O"), 2, 1))
ggplot(mfull_est, aes(reorder(term, term2), estimate, color = group, fill = group, shape = scale)) +
geom_hline(yintercept = 0, linetype = 2, color = "gray40", alpha = 0.5) +
geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0) +
geom_point(size = 2.5) +
scale_color_manual(values = pal2, name = "") +
scale_shape_manual(values = c(15, 19, 21, 17), name = "") +
scale_fill_manual(values = rep("white", 9), name = "") +
labs(y = "Estimate", x = "Standardized coefficient") +
scale_x_discrete(breaks = levels(mfull_est$term),
labels = c(expression(Herring[rec]),
expression(Herring[sub]),
expression(Sprat[rec]),
expression(Sprat[sub]),
expression(Cod[rec]),
expression(Cod[haul]),
expression(Flounder[rec]),
expression(Flounder[haul]),
expression(Saduria~entomon[rec]),
expression(Saduria~entomon[haul]),
expression(Depth[rec]),
expression(Depth[haul]),
expression(Oxygen[rec]),
expression(Oxygen[haul]),
expression(sigma[E]),
expression(sigma[O]),
expression(Temp[rec]),
expression(Temp[haul])
)) +
coord_flip() +
theme_plot() +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
plot.margin = unit(c(0.2, 0.2, 0.4, 0.2), "cm")) +
guides(color = "none", fill = "none", shape = guide_legend(ncol = 4))
ggsave("figures/effect_size.pdf", width = 17, height = 17, units = "cm")
# Just a test to see the labels were alright
ggplot(mfull_est, aes(reorder(term, estimate), estimate)) +
ggplot(mfull_est, aes(term, estimate)) +
geom_hline(yintercept = 0, linetype = 2, color = "red", alpha = 0.5) +
geom_point(size = 2) +
labs(x = "", y = "Standardized coefficient") +
coord_flip()
# Now save this for comparison in sensitivity analysis script
write.csv(mfull_est, "output/mfull_est.csv")
pred_map <- predict(mfull, newdata = pred_grid)
pal <- wes_palette("Zissou1", 21, type = "continuous")
plot_map_raster +
geom_raster(data = filter(pred_map, depth < 120 & depth > 10), aes(x = X*1000, y = Y*1000, fill = exp(est))) +
scale_fill_gradientn(colours = rev(pal), name = "Le Cren's condition factor") +
facet_wrap(~ year, ncol = 5) +
theme_facet_map() +
geom_sf(size = 0.1)
ggsave("figures/supp/condition/all_years_condition_map_covars.png", width = 6.5, height = 8.5, dpi = 600)
# And a smaller map for selected years
lab_size <- 1.7
plot_map_raster +
geom_raster(data = filter(pred_map, year %in% c(1994, 2001, 2008, 2018) & depth < 120 & depth > 10),
aes(x = X*1000, y = Y*1000, fill = exp(est))) +
scale_fill_gradientn(colours = rev(pal), name = "Le Cren's condition factor") +
facet_wrap(~ year, ncol = 2) +
theme(legend.key.height = unit(0.1, "cm"),
legend.key.width = unit(0.5, "cm"),
legend.position = c(0.9, .036),
legend.direction = "horizontal",
legend.background = element_rect(fill = NA),
legend.text = element_text(size = 6),
legend.title = element_text(size = 7.4)) +
geom_sf(size = 0.1) +
annotate("text", label = "Sweden", x = xmin2 + 0.25*xrange, y = ymin2 + 0.8*yrange, color = "black", size = lab_size) +
annotate("text", label = "Denmark", x = xmin2 + 0.029*xrange, y = ymin2 + 0.43*yrange, color = "black", size = lab_size, angle = 75) +
annotate("text", label = "Germany", x = xmin2 + 0.07*xrange, y = ymin2 + 0.025*yrange, color = "black", size = lab_size) +
annotate("text", label = "Poland", x = xmin2 + 0.55*xrange, y = ymin2 + 0.1*yrange, color = "black", size = lab_size) +
annotate("text", label = "Russia", x = xmin2 + 0.95*xrange, y = ymin2 + 0.2*yrange, color = "black", size = lab_size) +
annotate("text", label = "Lithuania", x = xmin2 + 1*xrange, y = ymin2 + 0.47*yrange, color = "black", size = lab_size, angle = 75) +
annotate("text", label = "Latvia", x = xmin2 + 0.99*xrange, y = ymin2 + 0.68*yrange, color = "black", size = lab_size, angle = 75)
ggsave("figures/condition_map_covars.pdf", width = 17, height = 17, units = "cm")
# Plot spatiotemporal random effect
plot_map_raster +
geom_raster(data = pred_map, aes(x = X * 1000, y = Y * 1000, fill = epsilon_st)) +
scale_fill_gradient2() +
facet_wrap(~ year, ncol = 5) +
theme_facet_map() +
geom_sf(size = 0.1)
ggsave("figures/supp/condition/epsilon_st_map.png", width = 6.5, height = 8.5, dpi = 600)
# Plot spatial random effect
plot_map_raster +
geom_raster(data = filter(pred_map, year == 1999), aes(x = X * 1000, y = Y * 1000, fill = omega_s)) +
scale_fill_gradient2() +
geom_sf(size = 0.1)
ggsave("figures/supp/condition/omega_s_map.png", width = 6.5, height = 6.5, dpi = 600)
# Fit a linear model to each prediction grid of the estimate over time
# https://community.rstudio.com/t/extract-slopes-by-group-broom-dplyr/2751/7
time_slopes_by_year <- pred_map %>%
mutate(id = paste(X, Y, sep = "_")) %>%
split(.$id) %>%
purrr::map(~lm(est ~ year, data = .x)) %>%
purrr::map_df(broom::tidy, .id = 'id') %>%
filter(term == 'year')
# Plot the slopes
time_slopes_by_year2 <- time_slopes_by_year %>%
separate(id, c("X", "Y"), sep = "_") %>%
mutate(X = as.numeric(X),
Y = as.numeric(Y))
plot_map_raster +
geom_raster(data = time_slopes_by_year2, aes(x = X * 1000, y = Y * 1000, fill = estimate)) +
scale_fill_gradient2(midpoint = 0, name = "Slope (condition factor~year)") +
geom_sf(size = 0.1)
<img src=“condition_model_cf_files/figure-html/calculate”spatial trend”-1.png” width=“672” style=“display: block; margin: auto;” />
ggsave("figures/supp/condition/spatial_trend_condition.png", width = 6.5, height = 6.5, dpi = 600)
# Prepare prediction data frame
nd_dep <- data.frame(depth_sc = seq(min(d$depth_sc), max(d$depth_sc), length.out = 100))
nd_dep <- nd_dep %>%
mutate(year = 2003L,
year_f = as.factor(2003),
ln_length_cm = 0,
biomass_her_sc = 0,
biomass_her_sd_sc = 0,
biomass_spr_sc = 0,
biomass_spr_sd_sc = 0,
density_cod_sc = 0,
density_cod_rec_sc = 0,
density_fle_sc = 0,
density_fle_rec_sc = 0,
density_saduria_sc = 0,
density_saduria_rec_sc = 0,
depth_rec_sc = 0,
oxy_sc = 0,
oxy_rec_sc = 0,
temp_sc = 0,
temp_rec_sc = 0
)
# Predict from full model
p_cond_dep <- predict(mfull, newdata = nd_dep, se_fit = TRUE, re_form = NA)
cond_dep <- ggplot(p_cond_dep, aes(depth_sc, exp(est),
ymin = exp(est - 1.96 * est_se),
ymax = exp(est + 1.96 * est_se))) +
geom_ribbon(alpha = 0.4) +
geom_line() +
theme(aspect.ratio = 1) +
labs(x = expression(Depth[haul]))
# Prepare prediction data frame
nd_oxy <- data.frame(oxy_sc = seq(min(d$oxy_sc), max(d$oxy_sc), length.out = 100))
nd_oxy <- nd_oxy %>%
mutate(year = 2003L,
year_f = as.factor(2003),
ln_length_cm = 0,
biomass_her_sc = 0,
biomass_her_sd_sc = 0,
biomass_spr_sc = 0,
biomass_spr_sd_sc = 0,
density_cod_sc = 0,
density_cod_rec_sc = 0,
density_fle_sc = 0,
density_fle_rec_sc = 0,
density_saduria_sc = 0,
density_saduria_rec_sc = 0,
depth_sc = 0,
depth_rec_sc = 0,
#oxy_sc = 0,
oxy_rec_sc = 0,
temp_sc = 0,
temp_rec_sc = 0
)
# Predict from full model
p_cond_oxy <- predict(mfull, newdata = nd_oxy, se_fit = TRUE, re_form = NA)
cond_oxy <- ggplot(p_cond_oxy, aes(oxy_sc, exp(est),
ymin = exp(est - 1.96 * est_se),
ymax = exp(est + 1.96 * est_se)))+
geom_ribbon(alpha = 0.4) +
geom_line() +
theme(aspect.ratio = 1) +
labs(x = expression(Oxygen[haul]))
# Prepare prediction data frame
nd_tem <- data.frame(temp_sc = seq(min(d$temp_sc), max(d$temp_sc), length.out = 100))
nd_tem <- nd_tem %>%
mutate(year = 2003L,
year_f = as.factor(2003),
ln_length_cm = 0,
biomass_her_sc = 0,
biomass_her_sd_sc = 0,
biomass_spr_sc = 0,
biomass_spr_sd_sc = 0,
density_cod_sc = 0,
density_cod_rec_sc = 0,
density_fle_sc = 0,
density_fle_rec_sc = 0,
density_saduria_sc = 0,
density_saduria_rec_sc = 0,
depth_sc = 0,
depth_rec_sc = 0,
oxy_sc = 0,
oxy_rec_sc = 0,
#temp_sc = 0,
temp_rec_sc = 0
)
# Predict from full model
p_cond_tem <- predict(mfull, newdata = nd_tem, se_fit = TRUE, re_form = NA)
cond_temp <- ggplot(p_cond_tem, aes(temp_sc, exp(est),
ymin = exp(est - 1.96 * est_se),
ymax = exp(est + 1.96 * est_se))) +
geom_ribbon(alpha = 0.4) +
geom_line() +
theme(aspect.ratio = 1) +
labs(x = expression(Temperature[haul]))
# Prepare prediction data frame
nd_cod <- data.frame(density_cod_sc = seq(min(d$density_cod_sc), max(d$density_cod_sc), length.out = 100))
nd_cod <- nd_cod %>%
mutate(year = 2003L,
year_f = as.factor(2003),
ln_length_cm = 0,
biomass_her_sc = 0,
biomass_her_sd_sc = 0,
biomass_spr_sc = 0,
biomass_spr_sd_sc = 0,
#density_cod_sc = 0,
density_cod_rec_sc = 0,
density_fle_sc = 0,
density_fle_rec_sc = 0,
density_saduria_sc = 0,
density_saduria_rec_sc = 0,
depth_sc = 0,
depth_rec_sc = 0,
oxy_sc = 0,
oxy_rec_sc = 0,
temp_sc = 0,
temp_rec_sc = 0
)
# Predict from full model
p_cond_cod <- predict(mfull, newdata = nd_cod, se_fit = TRUE, re_form = NA)
cond_cod <- ggplot(p_cond_cod, aes(density_cod_sc, exp(est),
ymin = exp(est - 1.96 * est_se),
ymax = exp(est + 1.96 * est_se))) +
geom_ribbon(alpha = 0.4) +
geom_line() +
theme(aspect.ratio = 1) +
labs(x = expression(Cod[haul]))
# Prepare prediction data frame
nd_spr <- data.frame(biomass_spr_sd_sc = seq(min(d$biomass_spr_sd_sc), max(d$biomass_spr_sd_sc), length.out = 100))
nd_spr <- nd_spr %>%
mutate(year = 2003L,
year_f = as.factor(2003),
ln_length_cm = 0,
biomass_her_sc = 0,
biomass_her_sd_sc = 0,
biomass_spr_sc = 0,
#biomass_spr_sd_sc = 0,
density_cod_sc = 0,
density_cod_rec_sc = 0,
density_fle_sc = 0,
density_fle_rec_sc = 0,
density_saduria_sc = 0,
density_saduria_rec_sc = 0,
depth_sc = 0,
depth_rec_sc = 0,
oxy_sc = 0,
oxy_rec_sc = 0,
temp_sc = 0,
temp_rec_sc = 0
)
# Predict from full model
p_cond_spr <- predict(mfull, newdata = nd_spr, se_fit = TRUE, re_form = NA)
cond_spr <- ggplot(p_cond_spr, aes(biomass_spr_sd_sc, exp(est),
ymin = exp(est - 1.96 * est_se),
ymax = exp(est + 1.96 * est_se))) +
geom_ribbon(alpha = 0.4) +
geom_line() +
#coord_cartesian(ylim = c(-4.65, -4.4)) +
theme(aspect.ratio = 1) +
labs(x = expression(Sprat[sub]))
# Prepare prediction data frame
nd_sad <- data.frame(density_saduria_rec_sc = seq(min(d$density_saduria_rec_sc),
max(d$density_saduria_rec_sc), length.out = 100))
nd_sad <- nd_sad %>%
mutate(year = 2003L,
year_f = as.factor(2003),
ln_length_cm = 0,
biomass_her_sc = 0,
biomass_her_sd_sc = 0,
biomass_spr_sc = 0,
biomass_spr_sd_sc = 0,
density_cod_sc = 0,
density_cod_rec_sc = 0,
density_fle_sc = 0,
density_fle_rec_sc = 0,
density_saduria_sc = 0,
#density_saduria_rec_sc = 0,
depth_sc = 0,
depth_rec_sc = 0,
oxy_sc = 0,
oxy_rec_sc = 0,
temp_sc = 0,
temp_rec_sc = 0
)
# Predict from full model
p_cond_sad <- predict(mfull, newdata = nd_sad, se_fit = TRUE, re_form = NA)
cond_sad <- ggplot(p_cond_sad, aes(density_saduria_rec_sc, exp(est),
ymin = exp(est - 1.96 * est_se),
ymax = exp(est + 1.96 * est_se))) +
geom_ribbon(alpha = 0.4) +
geom_line() +
#coord_cartesian(ylim = c(-4.75, -4.4)) +
theme(aspect.ratio = 1) +
labs(x = expression(Saduria~entomon[rec]))
# cond_dep + cond_oxy + cond_temp + cond_cod + cond_spr + cond_sad
# plot_annotation(tag_levels = 'A')
#
# ggsave("figures/supp/condition/conditional_effects.png", width = 6.5, height = 6.5, dpi = 600)
p_cond_dep2 <- p_cond_dep %>%
mutate(variable = "Depth (haul)") %>%
dplyr::select(est, est_se, depth_sc, variable) %>%
rename(value = depth_sc)
#> mutate: new variable 'variable' (character) with one unique value and 0% NA
#> rename: renamed one variable (value)
p_cond_oxy2 <- p_cond_oxy %>%
mutate(variable = "Oxygen (haul)") %>%
dplyr::select(est, est_se, oxy_sc, variable) %>%
rename(value = oxy_sc)
#> mutate: new variable 'variable' (character) with one unique value and 0% NA
#> rename: renamed one variable (value)
p_cond_tem2 <- p_cond_tem %>%
mutate(variable = "Temperature (haul)") %>%
dplyr::select(est, est_se, temp_sc, variable) %>%
rename(value = temp_sc)
#> mutate: new variable 'variable' (character) with one unique value and 0% NA
#> rename: renamed one variable (value)
p_cond_cod2 <- p_cond_cod %>%
mutate(variable = "Cod (haul)") %>%
dplyr::select(est, est_se, density_cod_sc, variable) %>%
rename(value = density_cod_sc)
#> mutate: new variable 'variable' (character) with one unique value and 0% NA
#> rename: renamed one variable (value)
p_cond_spr2 <- p_cond_spr %>%
mutate(variable = "Sprat (sub)") %>%
dplyr::select(est, est_se, biomass_spr_sd_sc, variable) %>%
rename(value = biomass_spr_sd_sc)
#> mutate: new variable 'variable' (character) with one unique value and 0% NA
#> rename: renamed one variable (value)
p_cond_sad2 <- p_cond_sad %>%
mutate(variable = "S. entomon (rec)") %>%
dplyr::select(est, est_se, density_saduria_rec_sc, variable) %>%
rename(value = density_saduria_rec_sc)
#> mutate: new variable 'variable' (character) with one unique value and 0% NA
#> rename: renamed one variable (value)
conds <- bind_rows(p_cond_dep2, p_cond_oxy2, p_cond_tem2,
p_cond_cod2, p_cond_spr2, p_cond_sad2)
ggplot(conds, aes(value, exp(est), ymin = exp(est - 1.96 * est_se),
ymax = exp(est + 1.96 * est_se))) +
geom_ribbon(alpha = 0.4) +
geom_line() +
facet_wrap(~variable, scales = "free_x") +
labs(y = "Le Cren's condition factor",
x = "Scaled value") +
theme_plot()
ggsave("figures/supp/condition/conditional_effects.png", width = 6.5, height = 6.5, dpi = 600)